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ABSTRACT 

We examine the Lagrangian divergence of the displacement field, arguably a more natural 
object than the density in a Lagrangian description of cosmological large-scale structure. This 
quantity, which we denote ip, quantifies the stretching and distortion of the initially homo- 
geneous lattice of dark-matter particles in the universe, ip encodes similar information as the 
density, but the correspondence has subtleties. It corresponds better to the log-density A than 
the overdensity 5. A Gaussian distribution in ip produces a distribution in A with slight skew- 
ness; in 6, we find that in many cases the skewness is further increased by 3. 

A local spherical-collapse-based (SC) fit found by Bernardeau gives a formula for ip's 
particle-by-particle behavior that works quite well, better than applying Lagrangian pertur- 
bation theory (LPT) at first or second (2LPT) order. In 2LPT, there is a roughly parabolic 
relation between initial and final ip that can give overdensities in deep voids, so low-redshift, 
high-resolution 2LPT realizations should be used with caution. The SC fit excels at predicting 
ip until streams cross; then, for particles forming haloes, ip plummets as in a waterfall to —3. 
This gives a new method for producing TV -particle realizations. Compared to LPT realiza- 
tions, such SC realizations give reduced stream-crossing, and better visual and 1 -point-PDF 
correspondence to the results of full gravity. LPT, on the other hand, predicts large-scale flows 
and the large-scale power-spectrum amplitude better, unless an empirical correction is added 
to the SC formula. 
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1 INTRODUCTION 

In the present, quite observationally successful theory of cosmol- 
ogy, the universe began with nearly uniform density everywhere. 
According to the theory of inflation, the small fluctuations in it be- 
gan as tiny quantum fluctuations, that 'inflated' to macroscopic size 
as the universe ballooned in its first instants. 

In an Eulerian description, the density and velocity fields at 
fixed comoving positions describe this process of structure forma- 
tion. In a Lagrangian description, on the other hand, the fundamen- 
tal object is the displacement field, a vector field measuring the co- 
moving distance particles have traveled from their initial positions. 
Fluctuations are not fundamentally in the density, but in the sepa- 
rations between particles. If the particles are arranged on a cubic 
lattice, as they often are in iV-body simulations, the density fluctu- 
ations are really deformations of this lattice. In underdense regions, 
the lattice stretches out; in overdense regions, it bunches together 
and forms structures. 

While the density is still relevant in a Lagrangian framework 
(as it sources gravity), the simplest scalar to construct from the dis- 
placement field (besides its magnitude, which is irrelevant for local 
physics) is its Lagrangian divergence. We denote this divergence 
as ip. It is the lowest-order invariant (with respect to rotations and 



translations) of the tidal tensor of the displacement field, and quan- 
tifies the angle-averaged stretching of the Lagrangian sheet in co- 
moving coordinates. Where a mass element becomes stretched out, 
tp increases, and its density decreases. For a potential displacement 
field (i.e. with zero curl), tp carries all of its information. 

Lagrangian dynamics have long been applied to cosmology, 
going back at least to Zel'dovich ( 1 970 ). It can be insightful to envi- 
sion the process of structure formation in terms of the dynamics of 
a Lagrangian 'sheet,' a viewpoint that for instance has been applied 
to classify the types of caustics (folds of this sheet) that can form 
(e.g. | Arnold et aL]|1982| |Arnold||2001). This viewpoi nt has got- 
ten some attention recentl y (e.g. |Shandarin et al.|2012| |Abel et aT] 
[20TT] |Falck et al.|2012b| |Neyrinck|2012^ . The 'sheer is initially 
flat in six-dimensional position-velocity phase space, with vanish- 
ing bulk velocity everywhere as the cosmic scale factor a — > 0. 
Seen in position space, gravity stretches out the sheet in underdense 
regions, and bunches it together in overdense regions. Assuming 
cold dark matter, the sheet never intersects itself in six-dimensional 
phase space, and instead folds up in rough analogy to origami. 

Several modifications to the original Zel'dovich approxima- 
tion (ZA) were proposed, including higher-order Lagrangian per- 
turbation theory (LPT). In LPT, streams of matter typically over- 
cross in high-density regions, failing to form the bound structures 
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that they would in full gravity. Several modifications of LPT (e.g. 
|Coles et al. 1993; Me lott et al.|1994| and references therein) have 
been proposed to solve this problem, another attempted solution for 
which also occurs in the present paper. The adhesion model, for in- 
stance IKofman & Shan darin 1988 , Ko fman et al.|1992|[Shan darin 
|2009| [Valageas & Bernardeau 201 1 ; Hiddi ng et al.|2012| > prevents 
stream crossing by introducing an effective viscosity, allowing the 
structure to be predicted in an elegant geometrical fashion. There 
are other ways of producing approximate particle realizations in a 
Lagrangian perturbation-theory approach, e.g. the PINOCCHIO al- 
gorithm ( Monaco et al. 2002]). 

One might argue that TV-body simulations have become com- 
putationally cheap enough that such approximate realizations have 
little use. For example, though, such techniques have proven quite 
useful in Bayesian initial-conditions reconstruction (Kitaura & An- 
|gulo||2012| |Jasche & Wandelt||2012| >, where full TV-body simula- 
tions would be far too slow. 

In this paper, we explore the behavior and properties of the 
Lagrangian spatial-stretching parameter ip. In Section|2] we review 
several approximations for ip in the literature and its relationship 
to the density field. In Section [3] we explore the relationship be- 
tween ip and density variables (the overdensity and the log-density) 
in a class of 'local Lagrangian' toy models, including a spherical 
collapse (SC) approximation that is most relevant to structure for- 
mation. In Section [4] we compare these approximations to results 
from an TV-body simulation, demonstrating the success of the SC 
approximation. In Section[5] we test a simple new way of produc- 
ing particle realizations using the SC approximation, and compare 
it to LPT approaches. 



2 APPROXIMATIONS FOR THE DISPLACEMENT 
DIVERGENCE 

There are several analytical approximations in the literature for the 
Lagrangian divergence of the displacement field. Where q is a La- 
grangian coordinate of a particle, we denote the displacement field 
as 'l'(q) = x(q) — q (with x the Eulerian position of a particle), 



and ?/j(q) 



*(q). Here (V 9 -) is the divergence operator 



in Lagrangian coordinates. Assuming that vE" is a potential field, 
V'(q) = Vq0(q), where <p is the displacement potential. All of 
the approximations used in this paper assume that ^ is a potential 
field, which implies that ip contains all of the information in ^Sf. 
In this section, we examine a few of these approximations. In full 
gravity, $ is not potential, i.e. it has a nonzero curl. However, as 
we show below, much of the large-scale clustering is captured with 
the potential-flow assumption. 



lagh & Jeong in prep), giving much-reduced 'transients' compared 
to ZA-produced initial conditions. 
At second order, 



T/>2LPT(q) 



+ D 2 V q 



2,(2) 



(2) 



where </>(q) is the total displacement potential, and D 2 is the 
second-order growth factor. D2(t) ~ — %D\ (r), the approxima- 
tion holding to better than 2.6 percent for 0.1 < fl m < 1 ( |Bouchet| 
|et al.|1993] >. These first- and second-order potentials are 



V^ (1) (q)=<Mq), 



v^ (2) (q) = E{^(q)^(q)-[^(q)] : 



(3) 
(4) 



We introduce here an isotropic, parabolic approximation to 
2LPT, i/>2LPT,parab, for '02LPT- Note that (suppressing the W su- 
perscripts) 



£< 



,-) = 



£< 



(5) 



Also note that in 3D, X^i 0,« 2 ' s bounded by |(V,0) 2 (in the 
isotropic case that <f> : u are equal for all i), and (V 2 ,^) 2 (in the 
case that <p t u — V 2 </> for some i, with all other components zero). 
Putting this together, with 1/6 (the isotropic case) ^ C2 ^ 1/2, 
and recalling that D 2 < 0, 



'(/'2LPT(q) = -D x 8 + D 2 
1 



**-£(«$)' 



> -D 1 S+^D 2 5 2 
1 2 

~ — Slin + -<5l in = ^2LPT,parab(q), 



(6) 

(7) 
(8) 



in the last line using the above approximation for D 2 . 

Fig.[T|shows the 2LPT mapping between ipi and ipf (V q • ^ 
at redshifts 49 and 0, respectively), as well as the log-density 
ln(l+<5), for particles in a ACDM 256 3 -particle, 200 ft -1 Mpc-box 
size TV-body simulation, analyzed and discussed further below. The 
particles were advanced using the 2LPT algorithm as described by 
|Scoccimarro| ( |1998] >. In this standard technique, also used in pro- 
ducing a ZA realization, first a Gaussian random field ipi is gener- 
ated, ipf is then estimated using an LPT approximation, and the fi- 
nal displacement field is generated in Fourier space with an inverse- 
divergence operator. 

For the top panel, the divergence was measured in Fourier 
space, the native technique in the algorithm that produced the par- 
ticle distribution, 



2.1 Lagrangian Perturbation Theory 



The Zel'dovich approximation (Zel'dovich 1970] ZA) is the first- 
order, linear approximation in Lagrangian perturbation theory 
(LPT). The ZA gives 



V*ta(q,*r) = -<5i in (q,r) 



DM 
£>i (to) 



<5iin(q, r ) 



(1) 



where 5n n is the overdensity linearly extrapolated with the linear 
growth factor D\, and to is some initial time. 

The second-order (2LPT) expression is more complicated, but 
still straightforward to implement. This slight added complexity 



seems worth the trouble for initial-conditions generation (Scocci 



marro 1998 , Crocce et al. 2006 , Tatekawa & Mizuno 2007 ; McCul 



(V • ij}) k = -tk ■ it> k 



(9) 



For the middle panel, and in the rest of this paper, was 
measured in real space by differencing Eulerian positions of par- 
ticles immediately before and after the particle at position q, in 
Lagrangian rows and columns of the initial lattice along the three 
Cartesian directions. 

There is a noticeable difference between the two methods of 
measuring ipf in the top panels. One reason for this is that the ef- 
fective resolution of the real-space ipf estimator is twice that of the 
Fourier-space estimator. 

Particularly using the Fourier-space estimator, the 2LPT pre- 
diction fails at high ip. Naively equating tpf and —5, this predicts 
strong overdensities in initial underdense regions ! This behavior is 
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Figure 1. Two-dimensional histograms showing the relation between the 



stretching parameter tp(<l) 



*(q), as well as the log-density, in the 



initial and final conditions, advancing particles according to second-order 
Lagrangian perturbation theory (2LPT). In the top panel, iji is measured 
in Fourier space; in the middle, it is measured by differencing particle po- 
sitions. The density in the bottom panel is estimated using a Voronoi tes- 
sellation, and is affected by multi-streaming. The color of each grid cell 
corresponds to the number of particles (out of 256 3 ) in that bin. In the top 
panels, dotted green lines show the linear relationship in the Zel'dovich ap- 
proximation, and the dashed red curve shows the isotropic approximation 

l/>2LPT,parab Of Eq.[8] 



tempered in the real-space-estimated ifif, but still there are many 
apparently overdense, initially underdense particles. 

This raises the question of whether artificial haloes might pop 
up in what should be voids using 2LPT. This is important since 
2LPT is sometimes used to generate low-redshift density distribu- 
tions, for example in modelling a sparsely sampled, large-volume 
survey, where only the clustering on large (e.g. baryon-acoustic- 
oscillation) scales needs to be accurate (e.g. Scoccimarro & Sheth| 



As one test of this issue, we show the density as well in Fig 
T] measured wi th a Voronoi tessellati on (e.g., |Schaap & van de 



Weygaert 2000; Neyrinck et al.|2005||van de Weygaert & Schaap 



2009[ >. For this Voronoi density estimate, each particle occupies a 
Voronoi cell, a locus of points closer to that particle than to any 
other particle. The overdensity <5vtfe = (V) /V — 1 at a particle 
is set according to the volume V of its cell. This density measure is 
mass-weighted, and thus in a sense Lagrangian, but only strictly La- 
grangian without multi-streaming, which does occur in this 2LPT 
realization. 

In the bottom panel of Fig. [T] at moderate to high densities, 
there is little correlation between tpi and 8. The over-shell-crossing 
in LPT, evident below in Fig.[l0] is one reason for this. At low den- 
sities, there are indeed a few overdense particles that have low 4>i. 
We find this also below in Fig. [12] after which we further discuss 
this issue. 

Previous authors ( |Buchert et~aT1|1994[ |Bouchet et al.||1995| 
Sahni & Shandarin 19961 have noted the failure of 2LPT in voids; 
they also found that going to third-order LPT (3LPT) does not 
improve agreement substantially. 3LPT comes at the expense of 
significantly greater complexity, and the addition of a nonzero 
curl component (as exists in full gravity, as well). Since we have 
adopted the approximation that the displacement field is potential 
in this paper, we stop our LPT analysis at second order. 

2.2 The Spherical-Collapse Approximation 



|Bernar deau ( 1994bjl gave a simple formula for the evolution of an 
average Lagrangian volume element, which gives a good fit to re- 
sults based on the spherical-collapse (SC) model. It is based on 
an Qm — > (and = 0) limit he (Bernardeau 1992) found to 
the nonlinear spherical-collapse evolution of density. A concise, in- 
structive derivation appears in (Bernardeau et al. 2002). Since the 
formula arises from a low-density limit, it is not surprising that it 
is quite accurate in voids in ACDM, so perhaps we should call it 
the spherical-expansion approximation. In this approximation, the 
mass element's volume, where Vb is the mean volume occupied by 
a particle (assuming equal masses), is 



V(t) = V 1 



T<5lir 



3/2 



(10) 



|2002{|Neyrinck & Sz apudi 2008l |Manera et al.|2012) T 



This approximation matches the behavior of 'rare events' (matter in 
deep voids) well. Fosalba & Gaztanaga ( 1998 1 discuss how this ap- 
proximation arises, to leading order, in a monopole approximation 
to Lagrangian perturbation theory. This expression has proven to be 
quite accurate, even outside of the 'rare-event' limit in which it was 
originally proposed, over a wide range of regimes and cosmologies 
l |Fosalba & Gaztanaga|1998) . |Scherrer & Gaztanaga] p00l) found 
that including a full parametric description of spherical-collapse 
dynamics further improved matters, however at the expense of ad- 
ditional complication and computation. 

To get an equation for the time evolution of ip out of this, we 
use a geometric isotropic-cube approximation to relate 8 and ip. 
The below derivation also essentially appears in Mohayaee et al. 
(2006). Assuming a Lagrangian mass element occupies a cube of 
side length -0/3 + 1 (giving V q ■ "3/ = tp), the volume of such a 
cube in units of the mean volume is 

V = 1/(1 + 5) = (1 + VV3) 3 . (11) 

Equating the RHS of Eq. \ \\ \ to the volume in Eq. [To] and 
employing the ZA V>iin = — <5un gives what we call the spherical- 
collapse (SC) approximation. 
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0.03 



V'sc = 3 



1 + gV>lh 



1/2 



- 1 



(12) 



1 /2 

6 

Ol(ro) 



where ^n n = j^j^jipo for some initial time To. This is one of a 
class of local approximations in which ipj at a given final redshift 
and position depends only on its linear value, ip\i n . In higher-order 
LPT, tpf is generally nonlocal, depending on derivatives of ipn n , as 
well. 



3 LOCAL LAGRANGIAN APPROXIMATIONS 

To explore some general relationships between the 'stretching pa- 
rameter' ip and density variables, we further explore them in a 
simple class of toy 'local Lagrangian' models introduced by |Pro-| 
|togeros &~S cherrer ( 1997 PS97). These models are parameterized 

by 1 < a < 3, 



5 a (i/j) = + 



(13) 



Here, ip is the actual tp of a volume element, not necessarily related 
to the linearly evolved ip\i n . 

It may help to think of a conceptually as the effective number 
of axes along which volume elements are expanding or contracting. 
The cubic-mass-element approximation in Eq. (11) has a = 3. 
However, confusingly, the a = 3 relationship was used in deriving 
the a = 3/2 SC approximation, a turns from 3 to 3/2 only when 
we add the spherical-collapse relationship of Eq. (10) , which relates 
('final') ip to the linearly evolved ijjwn- 

The a — 3/2 model is particularly useful, but a may take 
other effective values in other environments. So, we do not confine 
our attention exclusively to a = 3/2. 

In these models, density singularities arise at ip = —a, where 
the volume element has contracted to zero. In the SC approximation 
the critical density of a collapsed element is —lpn n = 1.5, close to 
the Einstein-de Sitter linear spherical-collapse density, 1.69. 

One way to quantify the non-linearity of the tp-5 relationship 
is in its Taylor-series coefficients. Eq. [13] expands to 



l + a 2 (l + a)(2 + a) 3 4 
-ip + _ W ~^y> V +0(ip ). (14) 



2a 



6a 2 



In this family of approximations, the log-density 
A a (ip) = ln(l + 5„) = -aln(l + ip/a) 



(15) 



has a much more linear relationship to ip than 5 does (recalling that 
a > 1): 



lot 



3q 2 



V> 3 + o(</> 4 ) 



(16) 



Curiously, the log-density is also closely related to the Eulerian 
divergence of the displacement field! Falck et al. 2012a), perhaps 
even more so than the Lagrangian divergence, as we investigate 
here 

In the next section, we will see that the distribution of A given 
a Gaussian-distributed ip is also significantly more Gaussian than 

5. 

At early epochs in the Zel'dovich approximation, a cubic- 
mass-element model with a = 3 describes the density distribution 
quite well. Fig. [2] shows the accuracy of this cubic-mass-element 
relationship is in a set of ZA-produced ACDM initial conditions 



0.02 



0.01 



0.00 




Figure 2. In green and black, two-dimensional histograms showing the re- 
lation between ip(q) = Vq ■ and <5, for 256 3 particles in a set of 
2 = 49 initial conditions produced using the Zel'dovich approximation 
(ZA). The deviation from tp = —8 follows the cubic-mass-element approx- 
imation <5 cu be for ip < 0; for ip > 0, the result is between <5 cu b c an d <5sc- 
The histogram is unnormalized, showing the number of particles in each 
bin. 



at redshift z = 49. This simulation, used below, has 256' 5 parti- 
cles, and box size 200 h^ 1 Mpc. Again, the density was estimated 
with a Voronoi method at each particle; here, it is a true Lagrangian 
density estimate, since the fluctuations are small enough that no 
multi-streaming has occurred. 

Although the nonlinearity in Fig.[2]is a bit accentuated by the 
stretched y-axis, it is still substantial. Putting a — 3 in Eqs. §14\ 
and \ 16) , the tp 2 coefficients in 5 and A are 2/3 (rather large; by 
far the highest among the approximations here explored) and 1/6. 



3.1 Density PDFs 

Analytical density PDFs easily emerge from such local Lagrangian 
approximations, which consist of simple transformations on initial 
distributions. Here we assume a Gaussian distribution in ip, but a 
non-Gaussian ip could also be transformed. 

Note that even without explicit initial non-Gaussianity at ar- 
bitrarily early times, a Gaussian distribution in ip results in a non- 
Gaussian S distribution. If a Gaussian S distribution is truly desired, 
one could start with an appropriately non-Gaussian ip distribution, 
though we do not explore this possibility here. 

We transform the distributions with the change-of-variables 
formula 



P(y) = P(x) \dx/dy\ , 



(17) 



where P(x)dx and P(y)dy give the PDFs of variables x and y. 

PS97 worked out the PDF of 5 for the above a-parameterized 
local-Lagrangian approximations, allowing shell crossing by using 
the absolute value of the volume element in Eq. (13) . Here, we take 
a slightly different approach, removing volume elements that have 
undergone shell crossing (i.e. with 1 + ip/a < 0) from consid- 
eration. This leads to a PDF that does not integrate entirely to 1, 
although its integral is bounded below by 1/2 for large cry,, and dif- 
fers negligibly from 1 for < 0.5. Assuming the fraction of such 
removed particles is small, i.e. for < 0.5, the PDF of ^Lag, the 
mass-weighted density distribution, is 



exp 



^ 2 {(l + 5)- 1 /«-l} 2 /(2^) 



(l + ^Lag)- 1 - 1 ^ 



(18) 



An Eulerian PDF for 5 can be obtained by multiplying the La- 
grangian PDF by a factor of V/(V), where V = 1/(1 + S), giving 
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Figure 3. Eulerian PDFs of 8 and A = ln(l + 5) from Eqs. {19} and (IT] , 

setting a = 3/2 and ctj, = 0.5. 
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5*0 from ZA 




ln(l+<5) 
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Figure 4. Skewness parameters 53, for both the overdensity 8 and the log- 
density A, letting — > 0, as a function of a, using Eqs. j!9| an d j21) . 
The curves are computed numerically, but match the relations in Eq. ]22[. 



exp 



P(foul) 



{(l + «5)- 1 /-_i} 2 /(2^) 



2na 2 



(19) 



The middle factor in the denominator is (V), which we found to 
have the form 1 + \ / 2(a)a 2 b for 1 < a < 3. The form for Vi 
gives the analytically calculable coefficient at V2(a = 1, 2, 3) = 
(1, 1/4, 1/3), and matches the numerically estimated coefficient at 
other q's, including the SC V 2 (3/2) = 1/6. 



3.2 Reduced non-Gaussianity in the log-density 



As is well-known in cosmology (e.g., |Coles & Jones|["l991| 
|Colombi|1994)|Neyrinck et al.|2009) , the PDF of the log-density, 
A = ln(l + S), is much more Gaussian than the PDF of S. One way 
to understand this is that P(S), unlike P (A), is tied down to zero at 
S = — 1. Because it is so easy to do in the a-parameterized model, 
here we give some explicit formulae for the PDF and skewness of 
A. 

A Gaussian tjj distribution transforms to 



exp 



P(A 



Lag ) 



-A/a 



2™$ 



(20) 



Figure 5. Volume-weighted skewnesses of A and <5, measured from 
Voronoi-cell volumes around particles in Zel'dovich realizations with 
power-law power spectra of indices n. The dotted curves are the asymp- 



values of S3 in the cubic-mass-element model. 



P(Abui) = 



exp 


■4, 


[e~ A ^ - l) 


2' 









(21) 



with the same (V) factor in the denominator as in Eq. 1 19 1. 

Fig. [5] shows Eulerian (volume-weighted) PDFs of S and A 
using the SC a = 3/2, with a^, = 0.5. Even at this modest a^, A 
is visibly more Gaussian than 5. 

The first-order non-Gaussianity statistic is the skewness S3 = 
(5 3 )/ (S 2 ) 2 , which has been worked out perturbatively in the mildly 
non-linear regime, in both Eulerian perturbation theory (EPT) and 
in the Zel'dovich approximation. Without any smoothing, to Eu- 
lerian second ('tree') order, S3 = 34/7 w 4.86 in an Einstein- 
de Sitter (EdS) universe (Peebles 1980), with small corrections 
in the ACDM case. In the Zel'dovich approximation, Sf cl = 4 
l |Bernardeau|1994a|[Fry & Scherrer|1994) , a bit lower. If the skew- 
ness measurement is done smoothing over equal-sized Eulerian 
cells, a term is added that depends on the (local) power-spectrum 
slope n c ff = dlna 2 (R)/dlnR, where R is the smoothing radius. 
With smoothing, 7 = — (n e g + 3) is added to S3. 

In several cases, we found that 53 in the limit of small fluctua- 
tions was reduced by 3 when measuring it from A instead of S. The 
simplest example is the exact lognormal distribution, for which, an- 
alytically, the skewness of the log-density S3 = 0, and S3 = 3, 
for any a a- A more general reduction of the skewness by 3 may 
only hold precisely in other cases in the limit a^, — > 0. 

Fig. H] shows our numerical (using Mathematica) estimate of 
5*3 and S3 as a function of a, letting a^, — > 0. The relations 



S3 = 3/q + 3, and S3 = 3/a 



(22) 



match the numerical solution, as well as our and PS97's analytical 
findings (at a — 1, 3/2, and 3). S3 (a = 3/2) = 5 is close to 
the full-gravity value from EPT mentioned above, 34/7 w 4.86, 
which is also the leading-order result in the full spherical-collapse 
model in the EdS case (Fosalba & Gaztanaga 1998, Bern ardeau] 
|et al. |2002) . In fact, S3 = 5 in the limit Q A ->■ of the full SC 
dynamics. This match to the a — 3/2 skewness is not surprising 
since the a — 3/2 model arises in the same limit. 
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It may be worth investigating tuning the a parameter further, 
for example to investigate a model with a = 21/13 « 1.62, which 
would exactly give the EdS EPT and SC skewness, and at the same 
time give a critical collapse ipn n = —1.62 for collapse, nearly the 
full nonlinear spherical-collapse value of ipu n « —1.686. 

As PS97 found, in the cubic-mass-element approximation as 
oy, — > 0, S3 = 4 as in the ZA, although S3 diverges from 4 dif- 
ferently than in the ZA as departs from 0. We numerically in- 
vestigated S3 in the ZA, measuring the volume-weighted particle- 
density skewness parameter S3 from several ZA realizations. The 
Zel'dovich-produced ACDM initial conditions of the simulation 
shown in subsequent sections have S3 = 4.01, and S3 4 = 1.07. 

Fig. [5] shows the volume-averaged skewness in 5 and A 
measured from the distribution of particle Voronoi densities in 
Zel'dovich simulations with power-law power spectra. The error 
bars are the dispersions among 3 realizations at each n. As n de- 
creases, large-scale over small-scale fluctuations dominate. S3 di- 
verges somewhat from 4 at high n. In the context of the a model, it 
makes intuitive sense that the isotropic, cubic-mass-element would 
be most valid for low n, where large-scale fluctuations dominate. 
As n increases, mass elements cease to expand or contract isotrop- 
ically, so the effective a decreases. 

Note that this measurement, although it is Eulerian (volume- 
weighted), does not include a smoothing of the type that would 
add a 7 term to S3, since there is no fixed cell size. Of course, the 
realizations have finite (mass) resolution; thus 'no smoothing' is 
not meant to imply infinite spatial and mass resolution. A 7 term 
from a fixed Eulerian cell size would cause S3 to depart from 4 in 
the opposite way than we observe when n is increased from -3. The 
used to generate each realization was held fixed at 0.02, and the 
power-law index n was varied from —1 to —2.5. Error bars show 
the dispersion from three different realizations analyzed at each n. 

From these numerical results, it appears that in the ZA, as well 
as in the a approximation, a log transform reduces S3 by 3. It would 
be interesting to show how widely this property holds, a question 
for later work. 



4 BEHAVIOR OF ip IN FULL GRAVITY 

Here we compare these theoretical estimates to what actually oc- 
curs in an TV-body simulation. The simulation has 256 3 parti- 
cles in a 200- h^ 1 Mpc box, run with a vanilla ACDM cosmology 
(fi m = 0.3, Q A = 0.7, cr 8 = 0.9, h = 0.73, n a = 1). The initial 
conditions were generated at redshift z — 49 using the Zel'dovich 
approximation, and run with the GADGET 2 ( Springel 2005 1 code. 
The spatial stretching parameter ip{n) = V 9 • vE'(q) is measured 
by numerically differencing neighboring particle positions. 

Fig.[6]shows the evolution of i/>(q) with redshift, showing 2D 
histograms of initial to final ip at different snapshots. At moderate 
and low densities, the initially straight line at z = 49 grows bent, 
and develops a scatter (accentuated somewhat by the color scale). 
At high densities, there is a critical value ip — —3, signifying col- 
lapse of the mass element. In an idealized case where a Lagrangian 
patch of particles contracts into a single point (a 'halo'), ip = —3. 
Here, the Lagrangian divergence of the particle-position field x(q) 
is zero; ip — V 9 ■ x — 3 = —3. 

The locality of the ipf-ipi relationship can be measured with 
the dispersion in these 2D histograms. At z — 0, as expected, the 
locality is high in void regions, and degrades at higher densities, 
but, considering the stretched color scale, the relationship is still 
rather tight. 
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Figure 6. In green, two-dimensional histograms showing the relation be- 
tween ipi at the initial redshift of 49, and ipf, measured at the redshifts 
listed. In a 'local Lagrangian' approximation, ipj is simply a function of 
rpi . With time, the dispersion between the two grows, but even at 2 = 0, 
the dispersion is rather low (considering the stretched color scale). Various 
local Lagrangian approximations are shown in magenta. The 2LPT curve, 
from Eq.[8] is a lower bound on V'2LPT- F° r a high-density mass element, 
ip f migrates downward until the element collapses, giving ip t = — 3, about 
which it oscillates afterward. The small white circles and lines about (0,0) 
show the magnitude of the shift in the approximation curves caused by en- 
forcing (ipf) = 0. 
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Figure 7. Trajectories in ip, as a function of cosmic scale factor a, of four sets of particles, 25 in each panel, drawn from the two-dimensional slice shown 
below of an iV-body simulation. From left to right: (1) A random set of particles; (2) particles within a Lagrangian distance of 2.3 h~ 1 Mpc of the maximum- 
ipi (lowest-initial-density) particle in the slice; (3) particles within a Lagrangian distance of 2.3 hT 1 Mpc of the minimum-»/>i (highest-initial-density) particle 
in the slice; (4) a random subset of particles with ipf < —3. Dotted curves show ip B c(ipi), using Eq. {12} , colored the same as each particle's actual trajectory. 
The dashed line indicates the ip = —3 collapse 'barrier.' 



As before, ip was measured by differencing the positions of 
Lagrangian-neighbor particles. The development of the ip = —3 
peak in ipf is sensitive to the method of measuring ip; it does 
not appear if the divergence is measured in Fourier space. Perhaps 
this arises from sharp edges being difficult to describe precisely in 
Fourier space. 

Once a particle crosses the ip = — 3 barrier, ip changes 
stochastically, but stays around —3, since Lagrangian neighbors 
stay nearby in a halo compared to the Lagrangian interparticle sep- 
aration (assuming somewhat low mass resolution). Thus a collaps- 
ing mass element's ip value evolves as though it were in a waterfall: 
it descends with time as the mass element's density increases, then 
hits a 'surface' at ip — —3, about which it then bobs around. 

For ipi < but ipf > —3 (at high densities before actual 
collapse), the best of the above approximations for ipf > —3 
seems to be V'2LPT,parab (Eq. ([8}), lying between the overpredict- 
ing Zel'dovich (Eq. {T}) and the underpredicting SC [Eq. \\2) ] pre- 
dictions. For ip ^ (at low densities), the SC prediction is best, 
again lying between the two rather poor alternatives. Curiously, 
these two approximations are both parabolic, one in ipi, and the 
other in ipf. 

We found empirically that another approximation shown, 

^ai fcxp = i3 1/2 (l-e- Bl/2 ^), (23) 

works well for both high and low ipi. However, we caution that 
to our knowledge it lacks theoretical motivation, and has strange 
behavior at very high ip (higher than plotted here), asymptoting to 

The curves do not precisely go through the origin, which ap- 
pears as a white dot. The offsets, shown by small white lines, en- 
sure that (ipse) = 0. For simplicity, we apply the same offset to all 
curves. The numerical value of this offset is similar for the various 
approximations, except for the ZA, for which ipf is always sym- 
metric about zero. This (ip) = condition ensures that there is no 
mean comoving expansion or contraction. 

The SC approximation, ijj sc , predicts a particular trajectory of 
ip with time, depending only on the local ipu n . To investigate how 
well this approximation holds with time, for particles in different 
environments, in Fig.lTlwe show trajectories in ip for particles in 



various classes: a random selection of particles; particles near the 
highest-initial-density and the lowest-initial-density particles; and 
particles with ipf < —3. At early epochs (a < 0.2), ip tracks ip sc 
well. Subsequently, particles participating in nonlinear structures 
can get seriously derailed (e.g., the rightmost panel). Still, there are 
many particles for which ip continues to track ip sc . In the deepest 
void, for example (second-left panel), the form of the expansion 
roughly tracks the SC prediction, but is skewed upward, perhaps 
due to the extremity of the void. 

Fig. [8] shows ipf, measured at z = for 256 2 particles oc- 
cupying a flat Lagrangian sheet from this simulation. Some similar 
figures appear in Mohayaee et al. (20061. The quantities are plot- 
ted in Lagrangian (initial-conditions) coordinates, with each pixel 
corresponding to a particle on the square lattice. Also plotted are 
the ORIGAMI {Falck et al.|20 12b l morphologies of the particles in 
the sheet, as well as the result of applying the spherical-collapse 
equation l | 1 _| i to the initial conditions. A particle's morphology in 
the ORIGAMI algorithm measures the number of Lagrangian axes 
along which other particles have crossed it in Eulerian space. A 
halo particle has been crossed by other particles along 3 orthogonal 
axes; filament, wall, and void particles have been crossed along 2, 
1, and orthogonal axes. In the bottom panels, if ipn n < —3/2, 
Eq. ^12\ has no solution, i.e. the mass element has collapsed; in 
this case, we set ip sc — —3. Contours indicate the boundaries of 
ORIGAMI halo regions. 

The color scheme suggests a topographical analogy, when 
working in Lagrangian coordinates: as time passes, ip departs from 
zero, in a way largely prescribed by its initial value. However, in 
overdense regions where it is decreasing, it is not allowed to plum- 
met arbitrarily; where collapses occur, 'lakes' form, where ip be- 
comes ~ — 3. 

In the upper-right panel, the blue 'lakes' of ipf correspond 
quite well to halo regions as identified by ORIGAMI. A simple halo- 
finder comes to mind, connecting particles on the Lagrangian lat- 
tice with ip under some threshold, approximately —3. We did try 
a simple implementation of this, but had difficulty finding a sim- 
ple threshold to characterize all haloes, since there are roughly as 
many halo particles with ip > —3 as ip < —3. Still, we suspect a 
halo-finder along these lines could be quite successful. 
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Figure 8. Quantities measured on a Lagrangian sheet of 256 2 particles from a 256 3 -particle simulation of box size 200 h~ x Mpc, run to redshift 0. Each 
pixel corresponds to a particle. Upper left: the ORIGAMI morphology of the particle (void, wall, filament and halo particles are colored black, blue, yellow 
and red). Upper right: the Lagrangian divergence of the displacement field, ipf- The 'lakes' are Lagrangian regions that have collapsed to form haloes. 
Lower left: a prediction from the initial-conditions 4>i > using the spherical-collapse formula j!2| . ipf is set to -3 if ip\in < —3/2. Lower right: the difference 
between upper left and lower right. The red contours mark the boundaries of haloes as identified by ORIGAMI. 



The bottom panels compare i/> sc and ipf. As in Fig. [6] the two 
match quite well in void regions, but in high-density regions, the 
correspondence is rougher. The rms difference between tp sc and ip f 
in this simulation is 1.31. The agreement in overdense regions can 
be improved by using tp sc for ip > 0, but V , 2LPT, P arab for ip < 0. 
This reduces the rms difference to 1.24. These are to be compared 
to the standard deviation of tpf (i.e. the rms if it is approximated 
with its mean, 0), which is 2.0. 



Another interesting feature of this plot is that ipf — ip sc (tpi ) of- 
ten plummets (i.e. becomes large and negative) on the Lagrangian 
outskirts of haloes. There are a couple of possible reasons for this. 
The particles have been dragged into the halo at late times, and may 
not be overdense initially. Also, as they have just fallen into haloes, 
their cubic Lagrangian volumes have likely just been swapped. 
Thus when a particle first collapses, t[> generally overshoots —3, 
as also shown in Fig. [7] 
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Figure 9. Voronoi-measured histograms of A = ln(l + 8), including all 
particles, and also including only ORIGAMI-identified void particles, which 
have undergone no stream-crossing. The theoretical PDFs are as in Eqs. 
\20\ and {21) , applying a normalization correction equal to the fraction 
of ORIGAMI void particles in the simulation. The top histogram is 'La- 
grangian' as in mass-weighted, with each particle contributing equally. The 
bottom, Eulerian histograms are also measured using the Voronoi tessella- 
tion; they are simply the 'Lagrangian' PDFs multiplied by V = 1/(1 + 5). 



Fig. [9] compares Voronoi-measured, mass-weighted PDFs to 
those assuming that each volume element evolves independently, 
i.e. the PDFs of Eqs. l |18|21} . Again, the 'Lagrangian' PDF is 
not truly Lagrangian, since the density estimate for a particle in- 
cludes all other particles, not just its Lagrangian neighbors. Thus 
stream-crossing boosts each particle's density according to the lo- 
cally overlapping number of streams, populating the high-density 
'shelf that poorly matches the approximation. 

At low densities, the approximations match the measurements 
well. PS97 also found this, albeit in simulations without as much 
structure. However, a normalization correction was necessary for 
the agreement in Fig. [9] A greater fraction of particles than pre- 
dicted by the SC approximation leave the SC tracks to populate the 
high-density tail, as in Fig. [7] 

If the SC approximation were precisely accurate in describing 
both the density evolution up to stream crossing, and the fraction 
of particles whose Lagrangian volumes have collapsed, no normal- 
ization correction would be necessary in Eq. \20\ , since as noted 
earlier, we did not divide by the integral over the PDF to assure a 
PDF integrating to unity. At this value of <j^, Eq. \20\ integrates 
to 0.71, which can be calculated with a simple Erf expression giv- 
ing the fraction of particles with (1 + fV'Un) < (the critical 
—ipiin = 1-5, intriguingly near the critical density for collapse, 
1.69). 

Instead, we found that a smaller factor, 0.34, gave a good fit to 
the low-density tail. This matches the fraction of void particles in 
this simulation snapshot as measured by the ORIGAMI (Fal ck et ah] 
|2012bj ) algorithm. An ORIGAMI void particle has not been crossed 
by any other particle over the course of the simulation. We also 
show a curve showing the PDF of only void particles. The shape 
of this curve does not quite match that from the SC-approximation 



PDF; that is, even higher-density void particles are scattered a bit 
to higher densities. 

The agreement between all curves looks much better in the 
bottom, Eulerian panel. This is because the PDF is simply the PDF 
in the upper panel multiplied by the volume factor V/ (V) . In terms 
of the x-coordinate A, this is simply an exponential damping, e~ A , 
bringing up the left side of the curve, and suppressing the right side. 

If the amplitude of this normalization correction can be esti- 
mated or calibrated accurately, Eqs. j20|21fr seem to provide a con- 
venient estimate for the nonlinear density PDF, if a PDF lacking 
the true, more-populated high-density tail is adequate. 



5 PARTICLE REALIZATIONS FROM THE SC 
APPROXIMATION 

Given that lp sc tracks the evolution of ip well, we explored how 
well it would work explicitly to advance tpi to ipf using Eq. 
{Til l, additionally setting = —3 for collapsed particles where 
ipiin < — 1.5 (the singularity in Eq. l|12}). This is a simple replace- 
ment for Rtpi-ipf relationship in a Zel'dovich code, entailing only a 
fast additional step (with a bit less computational effort than 2LPT). 
This may be about as well as one can do with a local prescription 
giving i/)/ as a function of V'lin, with no dependence on its deriva- 
tives. (For comparison to a 2LPT prescription, see Appendix D2 of 
|Scocci marro ( 1998 )). Here are all of the steps in our SC procedure: 

(i) Generate a Gaussian random field from a linear power spec- 
trum. This becomes ipn n . 

(ii) For ip\- m ^ —1.5, set i[>f — —3. For t/j/ > 1.5, set ipf — 
i'sc + C, where f/> sc is from Eq. (12) , and C is a constant ensuring 
that (tpf) — (easily measurable by summing up t/;/ with C = 0). 
C is typically small; for example, in Fig. [6] C is the length of the 
small white line attached to the white dot. 

(iii) Take the inverse divergence off/'/ (in Fourier space, invert- 
ing Eq. j9]l) to get the displacement field ty. 

Fig. [TO] shows the particle positions resulting from advanc- 
ing the initial (z = 49) conditions of this simulation to z = 0, 
compared to its actual z — particle positions. The SC approxi- 
mation gives a particle arrangement more visually similar to that 
using full A-body dynamics than using either of the LPT rela- 
tions, for example producing more concentrated 'haloes' where 
ip — —3. Still, they are not as pointlike as we had hoped, perhaps 
related to the difficulty of capture sharp edges in Fourier space. 
For example, measuring the divergence in Fourier space does not 
capture the ipf — —3 peak in Fig. [6] A completely real-space 
inverse-divergence algorithm might be more adept at producing 
tight haloes, but it is not obvious how such an algorithm would 
work. 

Fig. |1 1| shows mass-weighted 1 -point PDFs of each particle 
distribution. Here, the SC approximation gives the best approxima- 
tion to the true dynamics, especially for the lowest densities, where 
the agreement is also quite good in Fig. [6] There is a 'shelf of 
high-density particles; these correspond to particles in haloes. At 
higher resolution, this locus becomes a peak ( TFalck et al.|2012b) . 
The low-density peak in the 2LPT PDF is a sign of its inaccuracy 
at this rather high ^)n n dispersion, a(tf)n n ) — 2.7. 

Fig. [12] shows the same information in two-dimensional his- 
tograms, comparing the full-gravity (A-body) 8 to 5 in the approx- 
imately evolved realizations. Here, again ip sc performs best. In the 
low density region of the 2LPT scatter plot, there are in fact over- 
dense particles (reaching at maximum S ~ 8), which should be 
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Figure 10. Redshift-zero Eulerian locations of particles occupying a 256 2 sheet of a 256 3 -particle ACDM realization, projecting out the third dimension. 
Clockwise from upper left, particle positions are determined using: a full Af-body simulation; the Zel'dovich approximation; 2LPT; and the SC approximation 
j 1 2| By eye, the SC approximation gives the results closest to full gravity. 



extremely underdense in the final conditions. These middling over- 
densities are unlikely to produce spurious haloes detected in a 2LPT 
realization, but there remains some chance of that. 

What are we to conclude about the reliability of 2LPT at low 
redshifts for mock galaxy catalogs? The work here is hardly an ex- 
haustive study, as it considers just the single simulation analyzed 
here. But for this simulation at z = 0, the population of overdense 
particles that should be underdense starts to be a worry. This prob- 
lem would be even more severe if cr(V>im) were increased, popu- 
lating the high-f/>i branch of the i/>2LPT, P arab parabola. One way to 
increase a(tf)n n ) is by increasing the mass resolution (since fluctu- 



ations grow on small scales in a ACDM universe), so we recom- 
mend caution in using 2LPT realization at high resolution and low 
redshift. This is not surprising, of course; for high LPT accuracy, 
c(V'iin) should be < 1. Fortunately, to our knowledge, low-redshift 
uses of 2LPT have been at lower mass resolution than this, resulting 
in an appropriately low a(i/'iin)- 

While the SC approximation excels at predicting 1 -point 
statistics and a visually plausible particle distribution, unfortu- 
nately it seems to have deficiencies, as well. The SC approximation 
shifts the locations of nonlinear structures more than does the ZA. 
This is difficult to see in Fig. [TU] so we overplot the Af-body and 
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Figure 11. Mass-weighted histograms of Voronoi-estimated particle den- 
sities, for both the fully evolved initial conditions, and the three particle 
realizations shown in Fig. 1 10| The SC approximation performs best at low 
densities; all approximations fail to capture the second peak or shelf of high- 
density halo particles present in the full simulation. 2LPT even produces a 
PDF peak at low densities. 
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Figure 12. Two-dimensional histograms comparing particle densities 
evolved with full gravity (the rr-axis) to densities in the approximately 
evolved particle distributions. The dotted red lines show the ideal y = x 
locus. The SC-approximation-evolved particle distribution performs best. 
Note the turn-up at low densities in 2LPT, where some overdense particles 
are predicted among particles that should be in the deepest voids. 



approximate realizations in Fig.[T3] In 2LPT, structures have sim- 
ilar locations as in the ZA. We suspect that the discrepancy in the 
SC case is largely from voids that collapse in the full TV-body case. 
In the LPT approximations, overdense structures surrounding the 
doomed voids collapse as they should, but this is suppressed in the 
SC case. 

The differences in large-scale flows in the SC show up in root- 
mean-square errors in particle positions. The rms errors of particle 
positions compared to the full TV-body dynamics using the three ap- 
proximations are, respectively for Zel'dovich, 2LPT, and SC, 1.61, 
1.65, and 2.17 /i^Mpc. 

The LPT approaches are also more successful than the SC ap- 
proximation in predicting the low-redshift dark-matter power spec- 
trum amplitude on large scales, as shown in Fig. [14] The power 
spectrum of particles displaced according to the SC approximation 
gives a multiplicative bias on large scales of about 0.65 in this sim- 
ulation, although the SC power spectrum's shape is a bit closer to 
the shape of the full nonlinear power spectrum, turning down at 
smaller scales than do the LPT power spectra. 

To measure the power spectra in Fig. [14] particles were dis- 
placed according to each approximation, and then assigned to cells 
on a 256 mesh using Nearest Grid Point mass assignment. Power 
spectra were then measured from these meshes, correcting for shot 
noise. 

It is possible to fix this large-scale normalization issue by mul- 
tiplying ipn n by a factor in Eq. In Fig. [14] scaling i/j^ by an 
extra factor of 2.0 achieves this. The factor was found by iteratively 
changing the effective growth factor until the power spectra agreed 
on large scales. This 'corrected' approximation also improves the 
agreement between the SC and TV-body-evolved particles, bring- 
ing the rms error in particle positions down to 1.65 h~ Mpc, at the 
level of the LPT approximations. However, undesirably, the 'cor- 
rection' also reduces the SC power spectrum at small scales. 

In Fig. [15] we show the Fourier-space cross-correlation coef- 
ficient R(k) = Pg x s' /\/PsPs' between the z = density field 
in this TV-body simulation, and several other density fields. The 
'initial conditions' curve is essentially the propagat or (e.g. [ Crocce 
& Scoccimarro 2006). As pointed out recently by |Tassev & Zal- 
darriaga (2012), the cross-correlation between Zel'dovich and full 
gravitational dynamics is significantly higher at small k than the 
cross-correlation between Eulerian linear PT and the full dynam- 
ics. The agreement is even a bit better in 2LPT. The SC realization 
has poorer performance here, on the other hand, although with the 
bias correction, R(k) is the highest among the approximations at 
large k. 

Other ideas we tried in obtaining tp / locally from ipi were to 
use ip sc for ipi > 0, but V'2LPT,parab for tpi < 0, and also using 
VWfexp, based on the agreement in both cases to the data in Fig. [6] 
These results were quite similar to the simple ip sc approximation, 
though. 

It is quite likely that we could find some tweaking of these 
approximations that would empirically give the best results. How- 
ever, we did not explore this avenue exhaustively, since this best 
agreement could be limited to this single simulation. Still, it seems 
likely that further optimization of the ipi — >• tpf mapping would be 
fruitful. 



6 SUMMARY AND CONCLUSION 

In this paper, we examine tp, the Lagrangian divergence of the 
displacement field, arguably the most natural variable to quantify 
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Figure 13. The Zel'dovich and SC panels of Fig. |10| with the full A^-body results overplotted in red. While Zel'dovich gives artificially empty voids and 
fuzzier haloes, it gives somewhat more accurate large-scale flows than does SC. 




Figure 14. Matter power spectra in a 200- h^ 1 Mpc Af-body simulation 
at 2 = 0, compared to power spectra of particle distributions displaced 
according to various approximations. 



large-scale structure in a Lagrangian approach. The main results of 
the paper are as follows: 

• Even slight distortions of the initially uniform mesh of par- 
ticles, quantified by ip, produce a density distribution more log- 
normal than Gaussian. It seems that for a wide class of models, 
the skewness parameter S3 of the log-density field is reduced by 3 
compared to the skewness of the overdensity. 

• In 2LPT, the mapping from initial to final xp is roughly 
parabolic, allowing overdensities to form where there should be 
deep voids. This does not seem to be a significant worry for a mod- 
erately low-density redshift-zero 2LPT realization (< 1 particle per 
(ft -1 Mpc) 3 ), but caution is recommended at higher resolution. 

• The spherical-collapse-fit formula 1 ] 12f describes the evolution 



Figure 15. Fourier-space cross-correlation coefficients between the various 
approximately-evolved density fields and the particle distribution as evolved 
in the full TV-body simulation. The solid black line is essentially the non- 
linear propagator between the initial and final states; the Lagrangian cross- 
correlations are higher, indicating higher accuracy. 

of tjj from initial conditions better than first- or second-order La- 
grangian perturbation theory (LPT), up to halo formation. This also 
allows for an approximation to the 1 -point PDF of the density that 
works quite well for low-density, undisturbed (void) particles. 

• In LPT, ip gets arbitrarily small, indicating extreme stream 
crossing. In full gravity, however tj) gets stuck around —3, sig- 
nalling halo formation. This is the value it would have if Lagrangian 
regions contracted exactly to pointlike haloes. 

• This knowledge of how ip evolves allows for a new method 
to produce final-conditions particle positions, based on this SC ex- 
pression. Compared to LPT realizations, such SC realizations give 
reduced stream-crossing, and better visual and 1-point-PDF corre- 
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spondence to the results of full gravity. LPT realizations, on the 
other hand, give more accurate large-scale flows and large-scale 
power spectra, as well as improved cross-correlation to the den- 
sity field evolved with full gravity. An empirical correction may be 
added to the SC formula that seems to fix some of these issues, 
however. 

Our results suggest several possibilities for future work. We 
did not carefully investigate the new SC method of generating final- 
conditions particle positions with respect to redshift and resolu- 
tion. We suspect that the SC method could provide a good method 
of producing relatively low-redshift initial conditions for simula- 
tions, if such a thing is desired. It is true that the large-scale power- 
spectrum bias in this method is troubling, as are the shifts in large- 
scale flows, but these could be tied to ip = —3 collapses, and could 
be absent at high redshift without stream-crossing. 

The 'barrier' at ip — —3 could be useful for halo-finding in 
iV-body simulations. Unfortunately, it seems not straightforward 
to use this barrier to halo-find in a single snapshot of a simulation, 
but other possibilities exist. For instance, ip could be measured at 
each timestep; if a particle ever has ip ^ —3, it could be tagged as 
a halo particle. 

It is also quite interesting to consider ways of predicting where 
ip — —3 from the initial conditions. Such considerations may even 
allow provide analytical mass functions. Indeed, similar ideas have 
been proposed using LPT (e.g. Monaco et al. 20021. Another pos- 
sible approach may be to infer Lagrangian halo boundaries from 
V'scCV'iin) formula, as in the lower-left panel of Fig. [8] The true 
halo contours are often smoothed versions of these contours, and 
perhaps could be obtained by a combination of mathematical mor- 
phology techniques such as dilation and erosion (e.g. Serra 1983) 
in Lagrangian space, as can be useful in cleaning detected void 
boundaries ( |Platen et al.|2007fr . 

In conclusion, ip, a natural density-like variable in a La- 
grangian viewpoint, seems to be a rather useful quantity, with some 
extra information that is not in the density itself. It is fortunate that 
a simple formula gives ip's behavior in voids, where dark energy is 
most energetically dominant (if indeed it is a substance). To under- 
stand dark energy, understanding the stretching of the Lagrangian 
mesh in voids is likely particularly important. 
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